https://github.com/psotob91
Stata tiene varios comandos oficiales para realizar pruebas de hipótesis o estimar intervalos de confianza de manera clásica.
Tiene comandos para pruebas:
R base también tiene funciones para realizar inferencia estadística clásica.
Problema: No consistencia y dificulta de reporte reproducible.
Otros paquetes suplen este problema:
Provee marco simple e intuitivo compatible con tidyverse y pipelines.
Permite realizar varios tipos de pruebas:
Ingresa data.frame y retorna un data.frame
Múltiples funciones adicionales para remodelar, reordenar, manipular y visualizar los resultados de las pruebas.
Descargar e instalar paquete:
ttest y == numero
t_test(y ~ 1, mu = numero, detailed = TRUE)
Se usa función t_test() de paquete {rstatix}.
Permite obtener el valor p para una hipótesis dada y el intervalo de confianza para la media.
# A tibble: 1 × 12
estimate .y. group1 group2 n stati…¹ p df conf.…² conf.…³
* <dbl> <chr> <chr> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
1 60.5 weight 1 null model 50 -35.6 1.15e-36 49 58.3 62.8
# … with 2 more variables: method <chr>, alternative <chr>, and abbreviated
# variable names ¹statistic, ²conf.low, ³conf.high
El argumento mu = número es un valor de referencia contra el que deseamos compararnos.
El argumento detailed = TRUE permite obtener también los intervalos de confianza.
La función gt() es para mejorar visualización del resultado.
| estimate | .y. | group1 | group2 | n | statistic | p | df | conf.low | conf.high | method | alternative |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 60.524 | weight | 1 | null model | 50 | -35.57465 | 1.15e-36 | 49 | 58.29404 | 62.75396 | T-test | two.sided |
Agregaremos gt() a partir de ahora.
Interpretación:
La media estimada de glucosa fue de 95.48 mg/dL (IC95% 71.6 mg/dL a 119.4 mg/dL). Aunque este valor fue menor a 100 mg/dL (valor de referencia), no es posible concluir que la media poblacional es diferente del valor de refrencia 100 mg/dL (p = 0.627).
¿El nivel medio de estradiol es diferente de 50 UI en la población?
Para variable estradiol:
*No equivalente oficial aún.
sign_test(y ~ 1, mu = numero, detailed = TRUE)
La estructura es nombre_de_prueba_test
Prueba del Signo:
Permite comparar la mediana (ya no la media) contra un valor de referencia.
No permite estimar intervalos de confianza.
| estimate | .y. | group1 | group2 | n | statistic | p | df | conf.low | conf.high | method | alternative |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 58.8 | weight | 1 | null model | 50 | 50 | 1.78e-15 | 50 | 56 | 62.5 | Sign-Test | two.sided |
La mediana estimada de glucosa fue de 95.53 mg/dL. Aunque este valor fue menor a 100 mg/dL (valor de referencia), no es posible concluir que la media poblacional es diferente del valor de refrencia (p = 0.627).
¿El nivel mediano de estradiol es diferente de 50 UI en la población?
| estimate | .y. | group1 | group2 | n | statistic | p | df | conf.low | conf.high | method | alternative |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 103.95 | e2 | 1 | null model | 53 | 42 | 2.25e-05 | 53 | 75.78 | 122.97 | Sign-Test | two.sided |
La mediana estimada de estradiol fue de 103.95 UI. Este valor fue mayor a 50 UI (valor de referencia) (p < 0.001).
Se usa la prueba t de Student para datos pareados.
Como los datos están pareados 1:1, la prueba trabaja sobre la diferencia de los valores individuales convirtiendo el problema en una sola muestra.
Por lo tanto, los supuestos son los mismos que para la `prueba t de Student de un solo grupo.
Ejemplo: ¿Cuál es el efecto de la Vitamina C en el crecimiento de los dientes en cuyes experimentales?
len: Longitud de dientes al final del experimento.
supp: Suplemento de vitamina C (VC vs. OJ)
Error in file(file, "rt"): no se puede abrir la conexión
Error in mutate(., dif = VC - OJ): objeto 'df2' no encontrado
| estimate | .y. | group1 | group2 | n1 | n2 | statistic | p | df | conf.low | conf.high | method | alternative |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 3.7 | len | OJ | VC | 30 | 30 | 3.302585 | 0.00255 | 29 | 1.408659 | 5.991341 | T-test | two.sided |
Se tienen dos opciones: Prueba del signo vs. Prueba de rangos signados de Wilcoxon.
Se trabaja igual que para t-test pareada, con al diferencia de los valores de tal forma que se tiene solo un grupo.
Prueba del signo no requiere supuesto de simetría de distribución para comparar medianas.
Prueba de Wilcoxon requiere supuesto de simetría de distribución para comparar medianas.
Prueba del signo:
| estimate | .y. | group1 | group2 | n1 | n2 | statistic | p | df | conf.low | conf.high | method | alternative |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 4.85 | len | OJ | VC | 30 | 30 | 19 | 0.136 | 29 | -0.3 | 7.9 | Sign-Test | two.sided |
Prueba de Wilcoxon:
| estimate | .y. | group1 | group2 | n1 | n2 | statistic | p | conf.low | conf.high | method | alternative |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 4.000014 | len | OJ | VC | 30 | 30 | 575.5 | 0.0645 | -0.1000097 | 8.500017 | Wilcoxon | two.sided |
Supuestos de la Prueba del Signo
Aleatorización
Indenpendencia de observaciones
Escala de medición al menos ordinal
Supuestos de la Prueba de Wilcoxon
Aleatorización
Indenpendencia de observaciones
Escala de medición al menos ordinal
Distribución simétrica
prueba t de Student para grupos independientes.| estimate | estimate1 | estimate2 | .y. | group1 | group2 | n1 | n2 | statistic | p | df | conf.low | conf.high | method | alternative | p.adj | p.adj.signif |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| -2.2812500 | 59.21875 | 61.50000 | weight | Placebo | Dosis 1 | 16 | 16 | -0.7745439 | 0.445 | 30 | -8.296318 | 3.733818 | T-test | two.sided | 1 | ns |
| -1.5979167 | 59.21875 | 60.81667 | weight | Placebo | Dosis 2 | 16 | 18 | -0.7126840 | 0.481 | 32 | -6.164947 | 2.969114 | T-test | two.sided | 1 | ns |
| 0.6833333 | 61.50000 | 60.81667 | weight | Dosis 1 | Dosis 2 | 16 | 18 | 0.2249600 | 0.823 | 32 | -5.504008 | 6.870674 | T-test | two.sided | 1 | ns |
Supuestos de la Prueba t de Student para grupos independientes
Aleatorización
Indenpendencia de observaciones
Normalidad
Homogeneidad de varianzas (*)
| Name | Piped data |
| Number of rows | 53 |
| Number of columns | 14 |
| _______________________ | |
| Column type frequency: | |
| numeric | 1 |
| ________________________ | |
| Group variables | treat |
Variable type: numeric
| skim_variable | treat | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|---|
| weight | Placebo | 1 | 0.94 | 59.22 | 5.61 | 50.1 | 54.90 | 58.55 | 61.12 | 71.0 | ▅▇▆▃▃ |
| weight | Dosis 1 | 1 | 0.94 | 61.50 | 10.36 | 50.1 | 55.68 | 59.30 | 64.22 | 92.1 | ▇▇▁▁▁ |
| weight | Dosis 2 | 1 | 0.95 | 60.82 | 7.24 | 51.4 | 54.50 | 60.05 | 65.88 | 74.5 | ▇▂▃▃▃ |
| estimate | estimate1 | estimate2 | .y. | group1 | group2 | n1 | n2 | statistic | p | df | conf.low | conf.high | method | alternative | p.adj | p.adj.signif |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| -2.2812500 | 59.21875 | 61.50000 | weight | Placebo | Dosis 1 | 16 | 16 | -0.7745439 | 0.446 | 23.09487 | -8.372645 | 3.810145 | T-test | two.sided | 1 | ns |
| -1.5979167 | 59.21875 | 60.81667 | weight | Placebo | Dosis 2 | 16 | 18 | -0.7235647 | 0.475 | 31.45480 | -6.099330 | 2.903497 | T-test | two.sided | 1 | ns |
| 0.6833333 | 61.50000 | 60.81667 | weight | Dosis 1 | Dosis 2 | 16 | 18 | 0.2203118 | 0.827 | 26.44841 | -5.686973 | 7.053640 | T-test | two.sided | 1 | ns |
prueba de suma de rangos de Wilcoxon, también denominada U de Mann Withney.| estimate | .y. | group1 | group2 | n1 | n2 | statistic | p | conf.low | conf.high | method | alternative | p.adj | p.adj.signif |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| -32.14 | e2 | Placebo | Dosis 1 | 17 | 17 | 106 | 0.193 | -64.80 | 19.98 | Wilcoxon | two.sided | 0.579 | ns |
| -20.27 | e2 | Placebo | Dosis 2 | 17 | 19 | 125 | 0.257 | -68.78 | 26.52 | Wilcoxon | two.sided | 0.579 | ns |
| -0.40 | e2 | Dosis 1 | Dosis 2 | 17 | 19 | 160 | 0.975 | -41.24 | 33.54 | Wilcoxon | two.sided | 0.975 | ns |
Supuestos
Aleatorización
Indenpendencia de observaciones
Escala de medición al menos ordinal
Igual forma de distribución entre los grupos (solo diferencia en la posición)
Usaremos función anova_test() del paquete {rstatix}.
Esta función permite el uso de otros tipos de ANOVA: https://rpkgs.datanovia.com/rstatix/reference/anova_test.html
ANOVA Table (type II tests)
Effect SSn SSd DFn DFd F p p<.05 ges
1 treat 49.021 7171.016 2 100 0.342 0.711 0.007
# A tibble: 3 × 5
treat variable n mean sd
<chr> <fct> <dbl> <dbl> <dbl>
1 Dosis 1 weight 33 62.1 10.4
2 Dosis 2 weight 37 61.2 8.13
3 Placebo weight 33 60.4 6.40
Supuestos
Aleatorización
Indenpendencia de observaciones
Normalidad
+ Una mejor forma de evaluar la normalidad es usando los residuos.
Homogeneidad de varianzas (*)
Otra opción es evaluar los residuales. Esto lo veremos en regresión lineal.
Es preferible asumir que no hay homocedaticidad y realizar un ajuste robusto a heterogeneidad de varianzas.
ANOVA Table (type II tests)
Effect DFn DFd F p p<.05
1 treat 2 100 0.329 0.72
(*) Este supuesto puede flexibilizarse:
Usando algún método robusto de estimación de varianza.
Usando remuestreo (Prueba de permutación o Bootstrapping)
ANOVA Table (type II tests)
Effect SSn SSd DFn DFd F p p<.05 ges
1 treat 470.911 21640.35 2 103 1.121 0.33 0.021
# A tibble: 3 × 5
treat variable n mean sd
<chr> <fct> <dbl> <dbl> <dbl>
1 Dosis 1 lh 34 13.7 16.6
2 Dosis 2 lh 38 8.65 7.75
3 Placebo lh 34 11.6 17.7
Supuestos
Aleatorización
Indenpendencia de observaciones
Normalidad
+ Una mejor forma de evaluar la normalidad es usando los residuos.
Homogeneidad de varianzas (*)
Otra opción es evaluar los residuales. Esto lo veremos en regresión lineal.
A veces es preferible asumir que no hay homocedaticidad y realizar un ajuste robusto a heterogeneidad de varianzas.
ANOVA Table (type II tests)
Effect DFn DFd F p p<.05
1 treat 2 103 1.502 0.227
(*) Este supuesto puede flexibilizarse:
Usando algún método robusto de estimación de varianza.
Usando remuestreo (Prueba de permutación o Bootstrapping)
Alternativa no paramétrica de ANOVA one-way.
Si supuestos se cumplen, compara medianas.
Usaremos función kruska_test() del paquete {rstatix}.
Esta función permite el uso de otros tipos de ANOVA: https://rpkgs.datanovia.com/rstatix/reference/kruskal_test.html
# A tibble: 1 × 6
.y. n statistic df p method
* <chr> <int> <dbl> <int> <dbl> <chr>
1 weight 106 0.0347 2 0.983 Kruskal-Wallis
# A tibble: 3 × 5
treat variable n median iqr
<chr> <fct> <dbl> <dbl> <dbl>
1 Dosis 1 weight 33 60 9
2 Dosis 2 weight 37 61 12
3 Placebo weight 33 59 8.4
Supuestos
Aleatorización
Indenpendencia de observaciones
Misma distribución salvo por la mediana:
Esto significa que Kruskal Wallis también es perjudicado en cierto modo por la heterogeneidad de varianzas porque afecta el supuesto de igualdad de distribuciones
Podemos usar gráfico de densidades:
library(ggpubr)
datos %>%
ggdensity(x = "weight",
y = "..density..",
fill = "treat",
color = "treat",
add = "median")+ Podemos usar gráfico de violin: Combina cajas y densidad
Hay varios enfoques de comparaciones múltiples.
Algunas son pre-hoc y otras post-hoc.
Veremos solo un par para ejemplificar:
# A tibble: 3 × 9
term group1 group2 null.value estimate conf.low conf.high p.adj p.adj.sig…¹
* <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
1 treat Dosis 1 Dosis 2 0 -0.974 -5.80 3.85 0.881 ns
2 treat Dosis 1 Placebo 0 -1.72 -6.68 3.24 0.689 ns
3 treat Dosis 2 Placebo 0 -0.745 -5.57 4.08 0.928 ns
# … with abbreviated variable name ¹p.adj.signif
# A tibble: 3 × 8
.y. group1 group2 estimate conf.low conf.high p.adj p.adj.signif
* <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
1 weight Dosis 1 Dosis 2 -0.974 -6.39 4.44 0.903 ns
2 weight Dosis 1 Placebo -1.72 -6.85 3.42 0.701 ns
3 weight Dosis 2 Placebo -0.745 -4.91 3.43 0.904 ns
# A tibble: 3 × 9
term .y. group1 group2 df statistic p p.adj p.adj.signif
* <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
1 treat weight Dosis 1 Dosis 2 100 0.480 0.632 1 ns
2 treat weight Dosis 1 Placebo 100 0.824 0.412 1 ns
3 treat weight Dosis 2 Placebo 100 0.367 0.714 1 ns
# A tibble: 3 × 9
.y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
* <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
1 weight Dosis 1 Dosis 2 33 37 -0.0300 0.976 1 ns
2 weight Dosis 1 Placebo 33 33 -0.173 0.863 1 ns
3 weight Dosis 2 Placebo 37 33 -0.148 0.882 1 ns
Placebo Dosis 1 Dosis 2
Without couple 9 10 8
With couple 8 7 11
# A tibble: 1 × 6
n statistic p df method p.signif
* <int> <dbl> <dbl> <int> <chr> <chr>
1 53 1.04 0.594 2 Chi-square test ns
# A tibble: 6 × 9
Var1 Var2 observed prop row.prop col.prop expec…¹ resid std.r…²
<fct> <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Without couple Placebo 9 0.170 0.333 0.529 8.66 0.115 0.200
2 With couple Placebo 8 0.151 0.308 0.471 8.34 -0.118 -0.200
3 Without couple Dosis 1 10 0.189 0.370 0.588 8.66 0.455 0.789
4 With couple Dosis 1 7 0.132 0.269 0.412 8.34 -0.464 -0.789
5 Without couple Dosis 2 8 0.151 0.296 0.421 9.68 -0.540 -0.962
6 With couple Dosis 2 11 0.208 0.423 0.579 9.32 0.550 0.962
# … with abbreviated variable names ¹expected, ²std.resid
| treat | Total | p-value1 | |||
|---|---|---|---|---|---|
| Placebo | Dosis 1 | Dosis 2 | |||
| Marital status, recat | 0.6 | ||||
| Without couple | 9 (33%) | 10 (37%) | 8 (30%) | 27 (100%) | |
| With couple | 8 (31%) | 7 (27%) | 11 (42%) | 26 (100%) | |
| Total | 17 (32%) | 17 (32%) | 19 (36%) | 53 (100%) | |
| 1 Pearson's Chi-squared test | |||||
Placebo Dosis 1 Dosis 2
Without couple 9 10 8
With couple 8 7 11
datos2 %>%
tbl_cross(
row = married2,
col = treat,
percent = "row"
) %>%
add_p(
test = "fisher.test"
)| treat | Total | p-value1 | |||
|---|---|---|---|---|---|
| Placebo | Dosis 1 | Dosis 2 | |||
| Marital status, recat | 0.7 | ||||
| Without couple | 9 (33%) | 10 (37%) | 8 (30%) | 27 (100%) | |
| With couple | 8 (31%) | 7 (27%) | 11 (42%) | 26 (100%) | |
| Total | 17 (32%) | 17 (32%) | 19 (36%) | 53 (100%) | |
| 1 Fisher's exact test | |||||
library(gtsummary)
datos2 %>%
dplyr::select(treat, age, weight, e2, fsh, lh, race, married) %>%
tbl_summary(
by = treat
) %>%
add_p()| Characteristic | Placebo, N = 171 | Dosis 1, N = 171 | Dosis 2, N = 191 | p-value2 |
|---|---|---|---|---|
| Age, years | 34.0 (32.0, 37.0) | 31.0 (27.0, 38.0) | 35.0 (30.0, 36.5) | >0.9 |
| Weight, kg | 59 (55, 61) | 59 (56, 64) | 60 (54, 66) | >0.9 |
| Unknown | 1 | 1 | 1 | |
| Estradiol | 71 (40, 134) | 112 (91, 130) | 110 (76, 158) | 0.4 |
| Folicullo stimulant hormon | 3.66 (2.17, 5.49) | 4.00 (2.61, 4.83) | 4.37 (2.54, 5.30) | >0.9 |
| Luteinizant Hormon | 6 (3, 15) | 7 (4, 15) | 6 (4, 14) | >0.9 |
| Race | ||||
| Mestiza | 17 (100%) | 17 (100%) | 19 (100%) | |
| Marital status | 0.6 | |||
| Single | 8 (47%) | 9 (53%) | 8 (42%) | |
| Divorced | 0 (0%) | 1 (5.9%) | 0 (0%) | |
| Widowed | 1 (5.9%) | 0 (0%) | 0 (0%) | |
| Cohabiting | 0 (0%) | 1 (5.9%) | 0 (0%) | |
| Married | 8 (47%) | 6 (35%) | 11 (58%) | |
| 1 Median (IQR); n (%) | ||||
| 2 Kruskal-Wallis rank sum test; Fisher's exact test | ||||
No reporte valores p innecesarios, que no respondan sus preguntas de investigación pre-definidas.
Es una (muy difundida) mala práctica reportar tablas comparativas descriptivas con valores p.
Varias guías, comenzando por STROBE (para observacionales) y CONSORT (para ensayos clínicos), explícitamente recomiendan en contra de reportar estas tablas.
Solo haga una comparación descriptiva de los resultados de la tabla.
Reservese la inferencia para responder la pregunta principal:
Ahorra grados de libertad y previene grados de libertad "fantasma" (‘ghost degree of freedom’ según Frank Harrell)en contra de reportar valores p en las tablas descriptivas!!Fuente: Enlace aquí
mala práctica muy difundida: ¡Evitemosla en estudios observacionales!comando_modelo y x
funcion_modelo(y ~ x)
comando_modelo y x1 x2
funcion_modelo(y ~ x1 + x2, data = datos)
| Stata | R |
|---|---|
| y varnum varcat | y ~ varnum + varcat |
| y varnum, nocons | y ~ 0 + varnum |
| y i.varcat | y ~ varcat |
| y c.varnum1#c.varnum2 | y ~ varnum1:varnum2 |
| y c.varnum1##c.varnum2 | y ~ varnum1*varnum2 |
| y c.varnum1##i.varcat | y ~ varnum1*varcat |
regress y x1 x2 c.x3##i.x4 [weight = pesos]
modelo <- lm(y ~ x1 + x2 + x3*x4, data = datos, weights = pesos)
summary(modelo) # Opción de R base
modelo <- lm(y ~ x1 + x2 + x3*x4, data = datos, weights = pesos)
tidy(modelo) # Opción de R tidyverse con paquete broom
glm y x1 x2 [weight = pond], family(familia) link(enlace)
mod <- glm(y ~ x1 + x2,
family = familia(link = “enlace”),
data = datos, weights = pond)
summary(mod) # Opción de R base
tidy(mod) # Opción de R tidyverse con paquete broom
| Stata | Stata glm equivalente | R |
|---|---|---|
| regress | family(gaussian) link(identity) | family = gaussian(link = “identity”) |
| logit | family(binomial) link(logit) | family = binomial(link = “logit”) |
| poisson | family(poisson) link(log) | family = poisson(link = “log”) |
| nbreg | family(nbinomial ml) link(log) | Otra función: glm.nb() |
| NA | family(gamma) link(power -1) | family = Gamma(link = “inverse”) |
| NA | family(igaussian) link(power -2) | family = inverse.gaussian(link = “1/mu^2”) |
| NA | NA | family = quasi(link = “identity”, variance = “constant”) |
| NA | NA | family = quasibinomial(link = “logit”) |
| NA | NA | family = quasipoisson(link = “log”) |
ereturn list
names(mod) # Aunque lo mejor es explorar la lista guardada
predict cualq_nombre, tipo_residuo
cualq_nombre <- residuals(mod, type = “tipo_residuo”)
cualq_nombre <- augment(mod) # Genera lista de residuos y otros valores de interés
Usar un wrapper.
Stata tiene solo un empaquetado para regresión lineal con regress, pero aún no tiene empaquetados oficiales para glm u otras regresiones:
glm_diagnostic (Enlace aquí)no glm (cero inflados, beta y fraccionales, Cox y sobrevida paramétrica) y subirlo al repo oficial de Stata.R tiene infinidad de empaquetados, los más populares:
plot() de R baseautoplot() de {ggfortify}check_model() de {performance}Se puede usar broom para obtener resultados en data.frames y pipearlos dentro de paquetes para crear tablas personalizadas: flextable, gt, huxtable, kable.
collect.Se pueden usar empaquetados!
Realmente hay infinidad de empaquetados (wrappers) en R para tablas reproducibles.
Los más populares, estables y confiables son:
https://github.com/psotob91
percys1991@gmail.com
Curso: R para usuarios de Stata - Ed. 1 (2023)